## Create simulated data many times

simulateOutlier <- function(meth="ktimes", kfac = 2, stat="mean", covs="cont", ncovsexact=0, ncovscont = 2, mu =0, sd = 1, r=0, restr.list = NULL, n.tr=2, whenoutlier = NULL, howoutlier = "10timesmax", units=40, times=100, outfileDir=""){
	
	for (idx.times in 1:times){
	
    ## make bivariate normal data
    if(length(mu) < ncovscont){
      mmm <- rep(mu, ncovscont)
    }else{
      mmm <- mu
    }
	
    if(!(is.matrix(r))){
		cormat <- diag(1, ncovscont)
		cormat[row(cormat)!=col(cormat)] <- r
		SSS <- (sd^ncovscont)*cormat
		x <- mvrnorm(units, mmm, SSS)
	}else{
        xstan <- mvrnorm(units, rep(0, nrow(r)), r)
        x <- round(xstan*sd+mmm)
    }
	
    ## draw exacts from unif(cats)
    if(ncovsexact>0){
		for(idx.exdraw in ncovsexact:1){
			tmp.exdraw <- sample(restr.list[[idx.exdraw]], units, rep=TRUE)
			x <- cbind(tmp.exdraw, x)
			colnames(x)[1] <- cov.names[idx.exdraw]
		}
	}
		
    ## check draws given restrictions
    if(!(is.null(restr.list)) && (length(restr.list)>ncovsexact)){
		for(idx.unit in 1:units){
			for(idx.var in 1:ncol(x)){
				if(!(x[idx.unit, idx.var] %in% restr.list[[idx.var]])){
					x[idx.unit, idx.var] <- sample(restr.list[[idx.var]],1)
				}
			}
		}
	}

    if(covs=="contbim"){
		## make bivariate normal data
		mmm <- rep(mu, ncovscont)
		cormat <- diag(1, ncovscont)
		cormat[row(cormat)!=col(cormat)] <- r
		SSS <- (sd^ncovscont)*cormat
    consider.factor <- 50
		x <- mvrnorm(units*consider.factor, mmm, SSS)
		x <- x[order(x[,1]),]
		x <- x[c(1:(units/2), (units*consider.factor-units/2+1):(units*consider.factor)),]
	}

    cvnames <- paste("cont", 1:ncovscont, sep="")
	
    if(ncovsexact>0){
      exact.v <- cov.names[ncovsexact]
      exact.values <- x[1,1:ncovsexact]
    }else{
      exact.v <- exact.values <- NULL
    }
	
    if(!(is.null(whenoutlier))){
    if(howoutlier == "10timesmax"){
      for(idx.var in 1:ncovscont){
        wh.max <- which.max(abs(x[, idx.var]))
        x[whenoutlier, idx.var] <- 10*x[wh.max, idx.var]*sign(x[wh.max, idx.var])
      }
    }
  }
  
    seqblock1(id.vars = "id", id.vals =1, exact.vars = exact.v, exact.vals = exact.values, covar.vars = cvnames, covar.vals=x[1, (ncovsexact+1):ncol(x)], n.tr = n.tr, assg.prob.stat = stat, assg.prob.method = meth, assg.prob.kfac = kfac, file.name = paste(outfileDir, "/sim", idx.times,".RData", sep=""), verbose=F)
  
    for(idx.units in 2:units){
      seqblock2k(object= paste(outfileDir, "/sim", idx.times, ".RData", sep=""), id.vals=idx.units, covar.vals = x[idx.units, ], file.name = paste(outfileDir, "/sim", idx.times, ".RData", sep=""), verbose=F)
    }
	}
}
